Astronomy <fe Astrophysics manuscript no. 18346r 


© ESO 2012 


February 14, 2012 





Orbital characterization of the /3 Pictoris b giant planet* 

G. Chauvin^'^, A.-M. Lagrange^ H. Beust\ M. Bonnefoy^, A. Boccaletti^, D. Apai'', F. AUard^, 
D. Ehrenreich^, J. H. V. Girard^, D. Mouillet^, and D. Rouan^ 

^UJF-Grenoble 1 / CNRS-INSU, Institut de Planetologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, 
Grenoble, F-38041, France 

■^Max Planck Institute for Astronomy, Konigstuhl 17, D-69117 Heidelberg, Germany 
^LESIA, Observatoire de Paris Meudon, 5 pi. J. Janssen, 92195 Meudon, France 

* Department of Astronomy and Department of Planetary Sciences, The University of Arizona, 933 N Cherry Avenue, 
Tucson, AZ 85718 ^ Centre de Recherche Astronomique de Lyon, 46 allee d'ltalie, 69364 Lyon cedex 7, France 
^ European Southern Observatory, Casilla 19001, Santiago 19, Chile 



Received ; accepted 

ABSTRACT 

Context. In June 2010, we confirmed the existence of a giant planet in the disk of the young star /3 Pictoris located 
between 8 AU and 15 AU from the star. This young planet offers the rare opportunity to monitor a large fraction 
of the orbit using the imaging technique over a reasonably short timescale. It also offers the opportunity to study 
its atmospheric properties using spectroscopy and multi-band photometry, and possibly derive its dynamical mass by 
combining imaging with radial velocity data to set tight constraints on giant planet formation theories. 
Aims. We aim to measure the evolution of the planet's position relative to the star Pictoris to determine the planetary 
orbital properties. Our ultimate goal is to relate both the planetary orbital configuration and physical properties to 
either the disk structure or the cometary activity observed for decades in the j3 Pictoris system. 

Methods. Using the NAOS-CONICA adaptive-optics instrument (NACO) at the Very Large Telescope (VLT), we 
obtained repeated follow-up images of the /? Pictoris system in the Ks and L ' filters at four new epochs in 2010 and 
2011. Complementing these data with previous measurements, we conduct a homogeneous analysis, which covers more 
than eight yrs, to accurately monitor the /3 Pictoris b position relative to the star. 

Results. On the basis of the evolution of the planet's relative position with time, we derive the best-fit orbital solutions 
for our measurements. More reliable results are found with a Markov-chain Monte Carlo approach. The solutions favor a 
low-eccentricity orbit e ^ 0.17, with semi-major axis in the range 8-9 AU corresponding to orbital periods of 17-21 yrs. 
Our solutions favor a highly inclined solution with a peak around i = 88.5 ± 1.7°, and a longitude of ascending node 
tightly constrained at = —147.5 ± 1.5°. These results indicate that the orbital plane of the planet is likely to be 
above the midplane of the main disk, and compatible with the warp component of the disk being tilted between 3.5 deg 
and 4.0 deg. This suggests that the planet plays a key role in the origin of the inner warped-disk morphology of the 
/3 Pic disk. Finally, these orbital parameters are consistent with the hypothesis that the planet is responsible for the 
transit-like event observed in November 1981, and also linked to the cometary activity observed in the /3 Pic system. 
Conclusions. 

Key words. Techniques: adaptive optics, high angular resolution; Astrometry; Methods: observational, data analysis; 
Stars: planetary systems 



1. Introduction 

In the context of exoplanetary science, the direct imag- 
ing technique offers a unique observing window to explore 
the frequency and the properties of the extrasolar giant 
planet (EGP) population in large orbits (> 5 AU). The 
recent discoveries of massive planetary mass companions 
around stars and brown dwarfs (Chauvin et al. 2004, 2005; 
Lafreniere et al. 2008; Kalas et al. 2008; Marois et al. 2008; 
Lagrange et al. 2009; Marois et al. 2010) have shown that 
core-accretion alone cannot explain the formation of all im- 
aged giant planets, because the core-accretion timescales 
become too long and the disk surface density too low. 

Send offprint requests to: G. Chauvin 

* Based on observations collected at the European Southern 
Observatory, Chile (ESO programmes 072.C-0624, 384.C-0207, 
084.C-0739, 284.C-5057, 385.C-0132, 086.C-0341) 



Disk or core fragmentation are alternative mechanisms that 
could even operate in wide orbits, leading to different physi- 
cal and orbital distributions (Boley 2009; Dodson-Robinson 
et al. 2009; Vorobyov 2010). Additional mechanisms such 
as inward/outward migrations or planet-planet interactions 
might also modify the EGP orbital properties (Crida et al. 
2009). Consequently, there is a real need to characterize the 
orbital properties of recently imaged EGPs, to determine 
the planetary system architectures and dynamical stabili- 
ties, and to obtain additional constraints on their formation 
and evolution mechanisms. Follow-up astrometric studies 
have been conducted on the HR 8799 multiple exoplanetary 
system (Lafreniere et al. 2009; Marois et al. 2010; Bergfors 
et al. 2011; Soummer et al. 2011). Inclined orbital solutions 
for the b and d planets, and a 1:2:4 resonance for the planets 
c, d, and e, confirm stable configurations supported by dy- 
namical simulations (Reidcmeister et al. 2009; Fabrycky & 
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Fig. 1. Left: Astrometric positions of /? Pic b relative from A used in the present work to characterize the orbital parameters of 
the (3 Pic b planet. The predictions for the P Pic b in case of a background source are reported in gray from November, 10, 2003 
to March, 26 2011. Two linear fits, passing through the origin {Aa,A5) — (0,0), for the 2003 NE data points (at 34.25 deg) and 
the 2009, 2010 and 2011 SW measurements (at 211.15 deg) are indicated by the two blue dashed lines. A zoomed view of the most 
recent astrometric observations over 2010 and 2011 is presented on the right. 



Murray-Clay 2010). Owing to its even smaller semi- major 
axis, the planet (3 Pictoris b (hereafter (3 Pic b) offers a 
rare opportunity to constrain both the planetary orbital 
and physical properties (Lagrange et al. 2010, Currie et al. 
2011). This massive giant planet orbits at 8-15 AU around 
the young star l3 Pictoris (A5V, 19.44±0.05 pc, 121:4 Myr), 
which has been studied for almost three decades because of 
its emblematic debris disk (see Vidal-Madjar et al. 1998 for 
a review). So-called hot start theoretical models predict a 
mass of 9 ± 3 Mjup and T^s = 1700 ± 300 K (Quanz et al. 
2010, Bonnefoy et al. 2011). Characterizing the planetary 
properties offers a rare opportunity to directly investigate 
the planet-disk interactions, particularly the role played by 
the planet in the formation of the the inner warped compo- 
nent of the /3 Pic circumstellar disk (Mouillet et al. 1997; 
Augereau et al. 2001; Dawson et al. 2011). 



Since the recovery of /3 Pic b (i.e its re-discovery af- 
ter its passage behind the star), we have performed an as- 
trometric monitoring campaign, using NACO at the VLT. 
In Sect. 2, we describe observations aquired in the years 
2010 and 2011. In Sect. 3, we present the data analysis of 
these new data together with previous measurements in- 
cluding available astrometric calibrations (Lagrange et al. 
2009, 2010; Bonnefoy et al. 2011). Our main objective is to 
homogeneously analyze the astrometric position of /3 Pic 
b relative to the central star to accurately constrain its 
orbital properties. In Sect. 4, we report the most proba- 
ble solutions for the planet's orbit that best fits our as- 
trometric measurements, using two orbital fitting meth- 
ods, a least-square Levenberg-Marquardt algorithm and a 
Markov-chain Monte Carlo approach. Finally, in Sect 5, we 
discuss the consequences of our results in the context of 
previous astrometric studies, and their implications for a 
possible connection with the past transiting event observed 
in 1981, the disk structures, and the cometary activity of 
the P Pic system. 



2. Observations 

To monitor the (3 Pic b astrometry, we used the NACO 
high contrast adaptive optics (AO) imager of the VLT, 
equipped withthe NAOS AO system (Rousset et al. 2002), 
and the near-infrared imaging camera CONICA (Lenzen 
et al. 2002). The follow-up observations were obtained at 
five different epochs between September 2010 and March 
2011, using the angular differential imaging (ADI) mode 
of NACO. For accurate astrometry, two observing set-ups 
were used, the L ' filter with the L27 camera and the Kg 
filter with the S27 camera. The NACO detector cube mode 
was in addition used for frame selection. A classical dither- 
ing sequence was used with the repetition of five offset po- 
sitions to properly remove the sky contribution. In the end, 
the typical observing sequence represented a total of 200- 
250 cubes, i.e a total integration time of 35-50 min for an 
observing sequence of 1-1.5 hrs on target. The parallactic 
angles at the start and the end of our observations are re- 
ported in Table 1, together with the exposure time (DIT), 
the number of frame per cube (NDIT), and the number 
of cubes (Nexp) for each epoch. Typical exposure times 
of 0.15 s and 0.2 s were used in the Kg and L '-filters, re- 
spectively, to saturate the point spread function (PSF) core 
by a factor of 100 (a few pixels in radius) to improve the 
dynamics of our images. The observing sequences were ex- 
ecuted to optimize both the field rotation and the position 
of the secondary mirror diffraction spikes relative to the 
companion, except in March 26, 2011 when the data were 
obtained in service observing mode. Two sequences of non- 
saturated PSFs were acquired using a neutral density filter 
at the beginning and the end of each observing sequence 
to monitor the image quality. These data also served for 
the calibration of the relative photometric and astrometric 
measurements of /3 Pic b. During the different observing 
sequences obtained in visitor and service queue modes at 
ESO, the atmospheric conditions were stable, with seeing 
and coherence times of 0.6-1.0" and 3-6 ms , respectively. 
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Table 1. Obs log of new 2010-2011 VLT/NACO observations 



UT Date 


Filter/ Camera 


DIT 

(7 


NDIT 


cxp 


(dcK):(<lcg) 


28/09/10 


L7L27 
L'-ND;/L27 


0.2 


150 


199 


-60; -9 




0.2 


100 


20 


-62; -61 


16/11/10 




0.15 


100 


199 


-21; +14 




/s's-NDa/S27 


0.11 


100 


10 


-23; -21 


17/11/10 


L7L27 
L'-ND;/L27 


0.2 


100 


268 


-35; +32 




0.2 


100 


10 


-37; -35 


01/02/11 




0.15 


100 


223 


-6; +36 




A:,-ND,,/S27 


0.15 


100 


10 


-9; -6 


26/03/11 




0.15 


100 


249 


+38; +62 




/s:,-ND,/S27 


0.15 


100 


10 


-9; -6 



The same astrometric field was observed within a week of 
each follow-up observation of /3 Pic (see below). 

3. Data analysis 

For the present study, wc processed the data of the new 
observations of /3 Pic b obtained in September 28 2010, 
November 16 2010, November 17 2010, February 1st 2011, 
and March 26 2011. Previously archived data includ- 
ing available astrometric calibrations, obtained between 
November 2003 and April 2010, were also rc-proccsscd (sec 
Table 2 and Fig. 1 ). The most robust astrometric mea- 
surements (in terms of observing conditions and stability) 
were selected at each epoch. Data obtained on November 
16, 2010 and November 17, 2010 were both reduced to check 
the consistency of astrometric results obtained in the setups 
Z/'/L27 and Ks/S27 at a given epoch, as both setups were 
used for the astrometric analysis. The VLT/NACO M' data 
of Curric ct al. (2011) could not be considered owing to a 
lack of information about the NACO rotator offset position 
during the observation that could not be recovered from 
the ESO keywords, and that may affect the final absolute 
angular position. All data were homogeneously flat-fielded, 
cleaned from bad and hot pixels, and sky-subtracted. Sub- 
frames of 200 X 200 pixels were extracted to reduce the com- 
puting time of the data processing. The frames were recen- 
tered based on a registration of the central star position 
measured by a Moffat fitting of the non-saturated part of 
the stellar PSF wing, that was a threshold equal to at 1% 
of the detector linearity (i.e about 60% of the saturation 
limit). We also automatically decided to reject open-loop 
and poor- AO correction images, using a PSF encircled en- 
ergy criterion. 

For PSF-subtraction of the field-tracking data of 
November 11, 2003 and October 25, 2009, a reference star 
matching the /3 Pic observations in terms of parallactic an- 
gles was used to minimize the residuals in the final im- 
ages (see Lagrange et al. 2009). For ADI data, three differ- 
ent ADI-algorithms were applied to the PSF-subtraction, 
to verify the consistency of the solutions obtained using 
the classical ADI, smart ADI, and LOCI algorithms (see 
Lagrange ct al. 2010; Lafrenicre et al. 2007, for more de- 
tails on the algorithms). For smart ADI, we considered a 
separation criteria of 1.0 and 1.5 x FWHM, at the com- 
panion separation (ranging from 11 to 16 pixels), and ten 
images to compute each individual PSF. For LOCI, we con- 
sidered optimization regions of Na = 300 x FWHM, the 
radial-to-azimuthal width ratio ^ = 1, the radial width 



Ar = 2 X FWHM, and a separation criteria of 1.0 and 
1.5 X FWHM. Only smart ADI results were finally con- 
sidered owing the the ability of that method to reduce the 
companion self-subtraction, optimize the temporal PSF se- 
lection, and minimize the algorithm complexity when test- 
ing the photometric and astrometric systematics induced 
by the algorithm itself. We always found consistent results 
within the measurement uncertainty, which are detailed be- 
low, when using the classical ADI and LOCI algorithms. 

The main difficulty in deriving the planet's position rel- 
ative to the primary star was to accurately estimate the 
position of a central star with a saturated core, and the po- 
sition of the low signal-to-noise planet affected by the stellar 
residuals. To determine the central star position, as for the 
frames recentering, a Moffat fitting of the non-saturated 
part of the stellar PSF wing (with a similar threshold) was 
used. For the planet position and flux, we used a grid of 
5000 negative fake planets scanning a three-dimensional 
parameter space in (X,Y) positions (sampling of 0.02 pix- 
els) and flux (sampling of 10 ADU), injected one-by-one in 
the datacubes before PSF-subtraction. The datacubes were 
then reprocessed to derive the best-fit solution minimizing 
the residuals in the final subtracted-image, considering a 
region covering the companion ADI signature. The posi- 
tions of the companion relative to the primary star were 
finally transformed into sky coordinates using the 9i Ori C 
field observed with HST by McCaughrean & Stauffer (1994) 
(with the same set of stars TCC058, 057, 054, 034, and 026). 
For ADI data, the NACO rotator offset at the start of each 
ADI sequence was also calibrated and taken into account, 
using an astrometric binary observed in both field-tracking 
and pupil-tracking mode to estimate this offset. 

The results are given in Table 2 and reported in Fig. 1. 
The main contributors to the uncertainty in our astrometric 
measurements are listed below. Errors 2/ and 3/ were added 
quadratically, then linearly added to the systematic error 
in our measurement error 1/. The errors 4/ were neglected. 
We know consider the four types of error: 

1. The first contributor is the systematic error related to 
the determination of the (saturated) star center posi- 
tion, which is estimated from the fit to the PSF wings. 
We used non-saturated images to explore this effect 
given the saturation factor and the Moffat fitting thresh- 
old. If the PSF were centro-symmetric, the center esti- 
mate based on either the PSF wings or the PSF core 
would perfectly match. Tests on non-saturated data 
show that it does not. It induces a bias of 0.2-0.4 pix- 
els (i.e 6-12 mas). This asymmetry is variable such that 
this offset cannot be securely calibrated and subtracted. 
We note that, for consistency purposes, we re-analyzed 
the 2003 data with exactly the same method. Owing to 
the minor modification of the fitting pattern, this leads 
to slightly different values (of 0.5 pixels) from Lagrange 
et al (2009). 

2. A second source of error is the uncertainty in the com- 
panion position. Twelve negative fake planets of similar 
flux and separation at various position angles were then 
inserted to test the effect of both the stellar residuals in 
our measurements, and the measurement procedure it- 
self. The two non-saturated PSFs were used to take into 
account their influence on the flnal result. The typical 
error found was about 0.1-0.2 pixels (3-6 mas). 
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Table 2. NACO astrometric measurements of /3 Pic b relative to P Pic 



UT Date 


Mode 


Platoscalc 


True North" 


Rotator Offset' 


A a 


Ad 


Separation 


PA 




Ohs/Filtcr/Ohj 


(mas) 


(<l<'g) 




(mas) 


(mas) 


(mas) 


( 


10/11/03 


Fiold/L7L27 
Ficld/i /L27 
ADI/i7L27 


27.11 ±0.04 


0.29 ±0.07 


0. 


233 ± 22 


341 ± 22 


413 ± 22 


34.42 ± 3.52 


25/10/09 


27.11 ± 0.05 


-0.08 ± 0.10 


0. 


— 153 ± 14 


— 257 ± 14 


299 ± 14 


210.74 ± 2.89 


29/12/09 


27.10 ±0.04 


-0.06 ±0.08 


90.46 ± 0.10 


-163 ±9 


-260 ± 8 


306 ±9 


212.07 ± 1.71 


10/04/10 


ADI/i^,/S27 
ADI/L'/L27 


27.01 ± 0.04 


-0.26 ±0.09 


90.46 ±0.10 


-173 ± 7 


-300 ± 7 


346 ±7 


209.93 ±1.15 


28/09/10 


27.11 ± 0.04 


-0.36 ±0.11 


90.47 ±0.10 


-193 ±11 


-331 ± 11 


383 ± 11 


210.28 ±1.73 


16/11/10 


ADI/i(-,/S27 
ADI/L'/S27 


27.01 ± 0.05 


-0.25 ±0.07 


90.47 ±0.10 


-207 ± 8 


-326 ± 10 


387 ±8 


212.41 ± 1.35 


17/11/10 


27.10 ± 0.04 


-0.25 ±0.07 


90.46 ±0.10 


-209 ± 13 


-330 ± 14 


390 ± 13 


212.34 ±2.13 


01/02/11 


ADI/i^,/S27 


27.01 ±0.04 


-0.32 ±0.10 


90.46 ± 0.10 


-211 ±9 


-350 ± 10 


408 ±9 


211.13 ± 1.48 


26/03/11 


ADI/if,/S27 


27.01 ±0.04 


-0.35 ±0.10 


90.46 ± 0.10 


-214 ± 12 


-367 ± 14 


426 ± 13 


210.13 ± 1.81 



The orientation of true north is relative to the vertical of the detector, and is positive when lying to the east of the vertical. 
^ The NACO rotator offset position was calibrated and linked to the ESO keyword ADA.PUPILPOS according to the formulae 
ROTOFF = 179.44 — ADA.PUPILPOS, using various astrometric binaries observed at various epochs between November and 
December 2010. A conservative error of 0.10 deg was considered for this calibration. 



Table 3. Orbital solutions for /? Pic b: the best-fit reduced 
model obtained with the LSLM algorithm (top), and a typical 
"highly probable" orbit according to the MCMC fit (bottom). 
Note that we do not provide error bars here, as these are as- 
sumed to be described by the MCMC distribution. 



a(AU) P(yr) 


e 








tv (yr JD) 


Xr 


11.2 28.3 


0.16 


88.8 


-147.73 


4.0 


2013.3 


0.45 


8.8 19.6 


0.021 


88.5 


-148.24 


-115.0 


2006.3 


0.56 



3. We then considered the errors related to the platescale 
and true north orientation (see Table 2), as well as the 
uncertainty in the NACO rotator offset for ADI mea- 
surements (the error of 0.10 deg in the rotator angular 
accuracy). We note that an accurate absolute angular 
calibration of the NACO detector is difficult to achieve 
at mas precision. The two main limitations arc; the sig- 
nificant variation of the NACO detector true north with 
time, and the uncertainty associated with the calibra- 
tors themselves. Comparing data from different authors 
using different calibrators may then be risky in the end. 
A first solution is then to observe the same astrometric 
references at each epoch to ensure an accurate relative 
astrometry (as done here). A residual angular system- 
atic cannot be totally excluded but does similarly affect 
all measurements and the orbital solution (see discus- 
sion below). 

4. Finally, we neglected the errors related to the distortion 
correction, detector non-linearity, differential tilt jitter, 
Strehl ratio variation, or differential refraction which 
were estimated smaller than 1-2 mas; (see Fritz et al. 
2010). 

4. Orbital fit 

Our astrometric measurements reported in Table 2 show 
that the position angles are almost consistent with an edge- 
on orbital configuration, and are roughly consistent with 
the position angle of the /3 Pic circumstellar disk, which is 
reported to be 30.1 deg and 211.4 deg by Kalas & Jewitt 
(1995) for the NE and SW sides respectively, but with no 
error bars for these values). To derive the best- fit solution 



for our measurements, we considered the planet's inclina- 
tion as a free parameter. We assumed a Keplerian orbit 
described in a referential frame OXYZ^ where, as usual, 
the XOY plane corresponds to the plane of the sky and 
the Z-axis points towards the Earth. In this formalism, the 
astrometric position of the planet relative to the star is 
given by: 

X = = r (cos(a; ± w) cos f7 — sin(a; ± w) cosi sinJl) (1) 
y = Aa — r (cos(ti; ± w) sinfi ± sin(a; ± w) coszcosO) (2) 

where is the longitude of the ascending node (measured 
from North), a; is the argument of periastron, % is the incli- 
nation, V is the true anomaly, and r = a(l — e^)/(l±e cosu), 
where a represents for the semi-major axis and e the eccen- 
tricity. 

A Keplerian model was then fitted to our (Ai5, Aa) re- 
sults of Table 2, to derive the orbital period P (or equiva- 
lently the semi-major axis a, using the stellar mass Af* = 
1.75 Mq; Crifo et al. 1997), the e, i, 17, w, and the time 
for periastron passage t^. To constrain the /? Pic b's orbit, 
we used two complementary fitting methods: a least squares 
Levenberg-Marquardt (LSLM) algorithm (Press et al. 1992) 
to search for the model with the minimal reduced x^, and 
a more robust statistical approach using the Markov-chain 
Monte Carlo (MCMC) Bayesian analysis technique (Ford 
2005, 2006). The details of the MCMC analysis, which was 
adapted to this astrometric characterization, are reported 
in Appendix A. The best-fit LSLM solution found, and a 
typical example of a "highly probable orbit" according to 
the result of the MCMC study, are given in Table 3. These 
orbits are plotted in Fig. 2, in a plane containing the line 
of sight, as well as the positions of the planet at various ob- 
serving dates. The results of both orbital fitting methods 
are also shown for the five orbital elements (a, e, i, fi, w) in 
Fig. 3. In that figure, we also display the Xr distribution 
obtained for both analyses (considering a degree of freedom 
N — P = 12, where N is our number of measurements, here 
18, and P the number of parameters in our orbital model, 
here 6). Additional figures illustrating the goodness of both 
orbital fits are given in Appendix B. 

When comparing the results of the best-fit LSLM solu- 
tion with our MCMC distributions, our first striking result 
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Fig. 2. Plots ot the orbit of Table 3 with their orientation with 
respect to the line of sight. The larger orbit is the best LSLM 
model and the smaller one is an example of the highly prob- 
able orbit obtained with the MCMC approach. In each case, the 
dashed line shows the location of the periastron. The position of 
the planet at different observation epochs is shown as black dots 
along the orbit (projected error bars in the astrometric measure- 
ments are smaller than the symbol size) , and predictions for the 
upcoming years are shown in gray. 



is that they do not match. The best-fit LSLM solution has 
a relatively good Xr compared to the Xr distribution of the 
MCMC. However, this discrepancy seems to illustrate that 
our astrometric measurements are still too sparse to derive 
a deep LSLM Xr mininum, as the region of probable or- 
bits is probably fairly flat. Consequently, our confidence in 
the LSLM approach is low, and a robust determination of 
the orbital parameters of l3 Pic b is currently impossible 
with this fitting method. In contrast, we appear to have 
more success with the statistical MCMC approach, which 
produces the probable ranges of orbital elements of /? Pic 
b, enabling us to explore the significance of each orbital 
solution. 

The results of our MCMC analysis, reported in Fig. 3, 
are distributions of the orbital parameters that are far from 
gaussian, except for the inclination, and the longitude of 
ascending node. The semi-major axis distribution peaks at 
8.0-9.0 AU, and most eccentricities are given by e < 0.17. 
The inclination appears to be extremely concentrated close 
to 90°, but with a peak a.t i — 88.5 ± 1.7° (considering 
a confidence level at 1 a). This is indicative of a ~ 1.5° 
tilt with respect to a strict edge-on configuration. The lon- 
gitude of ascending node is also fairly well-constrained at 
n = -147.5±1.5° (f7 + 180 ^ 32.6±L5°, with a 1 cr error). 
Finally, the statistical distribution of the argument of peri- 
astron uj is more erratic, and does not appear to have clear 
solutions, which is not unsurprising owing to the low val- 
ues of the eccentricity distribution. Currie et al. (2011) per- 
formed a similar analysis, but on a smaller dataset that had 
not been homogeneously processed, thus with presumably 
larger systematics (different data processing, and possible 
lack of NACO rotator offset calibration). With the same 
dataset, we found similar results to Currie et al. (2011), us- 
ing our own MCMC code. The results of both studies, with 
four epochs (Table 1 from Currie et al. 2011) and eight 
epochs (Table 2, this work), are compared in Fig. 3. To a 
first order, they give similar results, although we note some 
differences. The peak of our semi-major axis distribution 



is shifted by about 1 AU to larger values. Our eccentricity 
distribution has a sharper cutoff that excludes a larger frac- 
tion of high eccentricity solutions, thanks to our most recent 
observations in 2010 and 2011. Our distributions of incli- 
nation, and longitude of ascending node are also shifted by 
respectively -1.0 and 1.5°. Finally, the distribution of Xr is 
better constrained in our case {Xr = 0.75±0.30, respectively 
to Xr = ^ 1-25), corroborating a more concentrated 
and robust set of orbital solutions. A peak value lower than 
unity however indicates that our uncertainty estimates are 
probably conservative. Further astrometric monitoring of 
/? Pic b, during the next quadrature, will be mandatory to 
improve the planetary orbital characterization. 

5. Discussion 

From our orbital fit analysis, four important outcomes arise: 
the semi-major axis of /? Pic b falls in the probable range 
of 8-9.0 AU, the eccentricity distribution is concentrated 
at e ^ 0.17, and the longitude of ascending node is fairly 
well-constrained at — 147.5 ± 1.5°, as is the inclination dis- 
tribution which peaks at 88.5±1.7°. The existence of a giant 
planet orbiting /3 Pic has been proposed by various studies 
during the past few decades (e.g. Freistetter et al. 2007). 
The main indirect indicators are (j) the inner warped com- 
ponent of the f3 Pic circumstellar disk, together with addi- 
tional asymmetries observed in the outer part (Mouillet et 
al. 1997; Kalas & Jewitt 1995); {ii) the photometric transit- 
like event observed in 1981 (Lecavelier des Etangs et al. 
1997); and (Hi) the cometary activity observed in the ab- 
sorption spectrum of (3 Pic (Ferlet et. al. 1987; Lagrange et 
al. 1996; Vidal-Madjar et al. 1998; Petterson & Tobin 1999). 
We discuss below how each of these observational findings 
may be related to both the existence and the orbital and 
physical properties of /3 Pic b. 

5.1. Disk - Planet configuration 

Dedicated scattered-light studies have accurately mapped 
the morphology of the (3 Pic disk (Kalas & Jewitt 1995; 
Heap et al. 2000; Golimowski et al. 2006). They revealed 
mainly a nearly edge-on disk composed of a main disk ob- 
served beyond 80 AU, and an inner warped component (at 
less than 80 AU), inclined by 2 — 5° with respect to the 
main disk position angle. The simulations of Mouillet et al. 
(1997), Augereau et al. (2001), and more recently Dawson 
et al. (2011) demonstrated that the presence of a planet 
orbiting the star at 10 AU, misaligned with the main disk, 
could actually form and sustain the /3 Pic inner warped 
disk. The observing challenge is then to test wether /3 Pic 
b might be this perturbating planet, that hence has with 
a significantly inclined orbit with respect to the main disk 
midplane. Currie et al. (2011) presented evidence of a mis- 
alignment between the planet and the inner warped disk 
of P Pic concluding that the planet was orbiting inside 
the main disk's orbital plane. We, however, do not confirm 
these results based on simultaneous measurements of the 
planet, and the main disk positions. This work, detailed in 
Lagrange et al. (2012), and using our i^s/S27 measurements 
of November 16, 2010 with a dedicated analysis for the disk 
orientation, shows that the current location of /3 Pic b is 
above the midplane of the main disk (which has a posi- 
tion angle of PAmain-disk = 209.0 ± 0.3 deg in the south 
west direction). This position is more compatible with the 
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Fig. 3. Results of the MCMC fit of the astrometric data of /3 Pic b: statistical distribution of the orbital elements a (top left), e 
(top middle), i (top right), Q (bottom left) and lo (bottom middle). We also show the distribution of Xr for the solutions obtained 
(bottom right). In each plot, the blue line show the results of our MCMC study. The dot-dashed line indicates the position of the 
best-fit LSLM Xr model obtained, and the -dashed line shows the position of an example of highly probable orbits according to 
our MCMC study see Table 3). The results of our MCMC study using the dataset of Currie et al. (2011) are indicated by a red 
line. 



warp component orientation tilted by 3.5-4.0 deg (i.e with 
P-A-warped-disk = 212.5 — 213.0 deg). Moreover, the current 
distribution of longitude of ascending node fi, which peaks 
at — 147.5 ± 1.5° (and would correspond to a position angle 
of 32.5 or 212.5°), can also be directly compared with the 
main disk and the warp orientations. This distribution cur- 
rently supports an orbital plane for (3 Pic b that does not 
coincide with the main disk midplane, but more probably 
with the warp component. The hypothesis that the /? Pic 
b planet is the one responsible for the inner warped mor- 
phology of the /3 Pic disk, without evoking any planetary 
inclination damping (an alternative scenario proposed by 
Dawson et al. (2011)), remains valid. We may therefore 
still conclude that the /? Pic b planet continues to excite 
the disk of planetisimals, forcing them to precess about its 
misaligned orbit. 



5.2. 1981 Transiting event 

Lecavelier des Etangs et al. (1995) reported significant 
photometric variations in November 1981 with a peculiar 
transit-like event on November 10, 1981. Lecavelier des 
Etangs et al. (1997) showed that a planet with 2-4 times 
the radius of Jupiter, orbiting at ~ 9AU at most could 
well be responsible for the photometric variations they re- 
ported. Lecavelier des Etangs & Vidal-Madjar (2009) in- 
vestigated this issue on the sole basis of the 2003 detection 
(Lagrange et al. 2009) of /? Pic b. They found that a transit 
of (3 Pic b in November 1981 could be compatible with a 
quadrature position in November 2003, assuming a semi- 
major axis in the range [7.6-8.7] AU. They also predicted 
that if the planet detected in 2003 on the NE side of the 
disk matched the one they predicted, it should reappear on 



the SW side in 2009 roughly where it was reobserved. We 
reinvestigate this prediction based on our MCMC orbital 
fit, considering both that we find a peak for the inclination 
distribution at 88.5 ± 1.7°, and the transit dates predicted 
in our set of orbital solutions. A tilt larger than 0.1 deg 
with respect to strictly edge-on configuration would sim- 
ply exclude the possibility that (3 Pic b is transiting along 
the line of sight. However, the sphere of influence of (3 Pic 
b with a Hill radius of about 1 AU (angular extension of 
about 7 deg at 8-9 AU) would still cross the line of sight 
(even considering all MCMC solutions), and therefore in- 
fluence the /? Pic photometry (if flUed with absorbing or 
scattering material). Assuming this hypothesis, we can still 
compare the date of the transit-like event of November 1981 
with the predicted transit dates of the sphere of influence 
of 13 Pic b between 1960 and 2030 obtained by our MCMC 
analysis (see Figure 4). The MCMC distribution of transit 
dates show that the parameters of the most recent (~ 1999) 
and next (~ 2018) transits are somewhat well-constrained. 
The peak is broader in ^ 1980. The suggested transit date 
of November 1981 falls to the right edge of that peak (al- 
though not at the center) and remains compatible with the 
current set of orbital properties obtained from our MCMC 
analysis. 



5.3. The j3 Pic cometary activity 

Transient redshifted spectral events have been regularly 
monitored in the absorption spectrum of /3 Pic (Ferlet et. 
al. 1987; Lagrange et al. 1996; Vidal-Madjar et al. 1998; 
Petterson & Tobin 1999), and attributed to the sublimation 
of numerous star-grazing planetesimals crossing the line of 
sight, which are referred to the falling evaporating bod- 
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Fig. 4. MCMC fit to tiie transit dates of /3 Pic b in front of 
the line of sight. The plotting conventions are the same as in 
Fig. 3. In addition, the date of the transiting event observed 
by Lecavelier des Etangs et al. (1997) is indicated by a thick 
vertical line. 

ies (FEBs) phenomenon. Their origin has been tentatively 
related to mean-motion resonances with a Jovian planet or- 
biting the star (Beust & Morbidelh 1996, 2000; Thebault 
& Beust 2001; Beust & Valiron 2007). Several constraints 
have been deduced from dynamical studies of the FEBs 
scenario, suggesting that: 

1. The planet responsible for the phenomenon is mas- 
sive enough (~ jovian) to allow numerous bodies to be 
trapped in mean-motion resonances. 

2. Its orbit is slightly eccentric (e ^ 0.05-0.1) to allow 
bodies trapped in the resonances to see their eccentricity 
pumped up (Beust & Morbidelh 1996, 2000). 

3. The longitude of periastron zu of the planet with respect 

to the line of sight is 70° ± 20° (Thebauh & Beust 

2001), to enable the statistics of the Doppler velocities 
of the FEB spectral signatures to match the observed 
ones which are largely strongly biased towards redshifts. 

4. Finally, the planet location is no further away than 
20 AU, otherwise the FEBs would hardly be able to get 
into the dust sublimation zone. 

The /3 Pic b planet obviously has orbital and physi- 
cal properties compatible with constraints 1/ and 4/. The 
situation is less straighforward for the constraints 2/ and 
3/. Eccentricities higher than ~ 0.05-0.1 are indeed fully 
compatible with our fit, but circular orbits cannot not ex- 
cluded. Finally, the longitude of perisatron w measured 
from the line of sight is related to the argument of peri- 
astron uj in our fit. The parameter uj is measured from the 
XOY plane of our reference frame, i.e., the plane of the 
sky. Assuming an edge-on orientation of the disk, we then 
have ui = zu + Tr/2. Thus, m ~ —70° ± 20° means that 
oj ~ 20° ± 20°. Unfortunately, our constraint on oj remains 
too low to state whether condition 3/ is fulfilled or not, 
partly due to our upper limit on the eccentricity distribu- 
tion at e ^ 0.17). Further measurements are needed to yield 
conclusive results. 

6. Conclusion 

We have reported the results of our analysis of follow-up 
observations of the astrometric positions of f3 Pic b rel- 
ative to /3 Pic. Together with previously archived data 
including available astrometric calibrations, we have ho- 
mogeneously reprocessed all relevant astrometric measure- 
ments of /3 Pic b at nine different observing epochs. We 



have taken into account the various contributors to the un- 
certainty in our astrometric analysis, including the deter- 
mination of the star position in the low-flux part of the 
saturated PSF, the influence of the PSF residuals on the 
companion position, and the errors related to the detec- 
tor calibration, and the NACO Nasmyth-rotator position. 
We then used orbital-fltting techniques to derive the most 
probable orbital solutions for the /3 Pic b planet: a least 
squares Levenberg-Marquardt algorithm and a Markov- 
chain Monte Carlo Bayesian analysis. As our measurements 
do not cover the complete planetary orbit, and are biased 
towards the most recent epochs since the planet recovery, 
the Markov-Chain Monte Carlo approach provides more ro- 
bust and reliable orbital solutions for /3 Picb. The most 
probable ones favor a low-eccentricity orbit e ^ 0.17, with 
a semi-major axis of 8-9 AU corresponding to orbital peri- 
ods of 17-21yrs, an inclination that peaks at 88.5 ± 1.7°, 
and a longitude of ascending node fairly well-constrained 
at —147.5 ± 1.5°. Our results support the previous astro- 
metric studies of Lagrange et al. (2010) and Currie et al. 
(2011), although our present study have indeed provided a 
more robust set of orbital solutions. Our conclusions sup- 
port the idea that the planet is not in the midplane of the 
main disk. The planet's position is more compatible with 
the warped component orientation, thus corroborating the 
idea that /3 Pic b is responsible for the inner warped disk 
morphology. The current orbital solution of /3 Pic b is still 
consistent with the sphere of influence of the planet being 
responsible for the 1981 transiting event, although a devi- 
ation of more than 0.1 deg from an edge-on configuration 
would exclude a planetary transit. Finally, the planet's pre- 
dicted mass, eccentricity, semi-major axis, and longitude of 
periastron also imply that /3 Pic b is likely to be the ori- 
gin of the cometary activity observed in the /3 Pic system. 
Further deep imaging characterization should help us to 
more tightly constrain the orbital parameter space, once the 
planet has passed the quadrature (most probably in 2013). 
Further spectroscopic or multi-photometric observations 
should also help us to determine the underlying physics of 
this giant planet in the framework of current planetary at- 
mosphere studies (Bonnefoy et al. 2010; Janson et al. 2010; 
Skemer et al. 2011, Madhusudhan et al. 2011; Barman et 
al. 2011a, 2011b). Therefore, much is to be expected from 
future extreme- AO instruments SCExAO/Subaru (Guyon 
2010), GPI/Gemini (Macintosh et al. 2008), SPHERE/VLT 
(Beuzit et al. 2008) in the coming years. 
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Appendix A: Convention and Markov Chain Monte 
Carlo adapted for Astrometry 

The astrometric position of the planet relative to the star 
is described by Eqs. (1) and (2). Fo any orbital solution, 51 
and uj changed to + tt, w + tt, u + tt respectively yields 
the same astrometric data. In the context of f2, the differ- 
ence between O and + vr actually resides in the z-motion 
(along the line of sight). If we consider a nearly edge-on 
orbit, we expect the longitude of the ascending node 51 to 
match the PA of the astrometric data (here, Vt ~ 34° or 
Q. ~ 211 = —149° if we consider an edge-on orbit for /3 
Pic b). By convention, 51 must thus corresponds to the PA 
when the z-coordinate grows, i.e., when the planet is mov- 
ing towards the observer. In the context of /3 Pic, the rota- 
tion sense of the gaseous disk was determined by Olofsson 
et al. (2001): the NE branch is receding from us while the 
SW branch is approaching. If we assume that the planet 
is moving in the same sense as the disk, then j3 Pic b was 
receding in the 2003 observations and is approching now, 
and it passed behind the star in between. This means that 
the ascending node is located towards the SW branch of the 
orbit, or that 51 ~ —149°. We will use this property to dis- 
tinguish between solutions that yield the same astrometric 
data. Note that this determination occurs in any case after 
the fitting procedure. In both approaches, we fit w -I- 51 and 
w — 51. Eqs. (1) and (2) can be rewritten unambiguously as 
a function of a; -I- f2 and w — 51). 

In the context of the Markov Chain Monte Carlo, let 
us briefly recall the detail of that technique that we here 
adapted for j3 Pic b. Let us call x the model parameters 
vector we want to constrain and d the data vector. We 
want to determine the posterior distribution p{x\d), i.e., 
the probability function of the parameter vector x given 
the data vector d and its error vector. This requires the 
knowledge of a prior probability function pq{x) for x. A 
Markov chain is a sequence of successive set of trial values 
a^i (i > 1) for a;. The Metropolis- Hastings (M-H) algorithm 
(Ford 2005) is used to derive aj^+i from Xi via a transition 
probability. After convergence, the equilibrium distribution 
of the chain equals the posterior distribution p{x\d). Ford 
(2006) suggests several assumptions for MCMC adapted to 
the search of exoplanets by radial velocity, depending on the 
kind of orbit we are looking for. We adapt here his recom- 
mendations to our astrometric fit. Following Ford (2006), 
we assume the prior distribution pq{x) to be uniform in 
X = (logP, e, cosi, H. + Lu,Lu — n, tp). However, the work is 
done on the parameter vector u{x). 

Our main motivation in using m as a variable instead 
of X was to improve the convergence of the Markov chains, 
as suggested by Ford (2006). A good way to achieve this is 
to avoid singularities. For instance, using e.cos{uj + 51) and 
e.sin{uj + instead of uj and 51 removes the non-regularity 
at e = (51 is not defined); using a; -I- 51 causes these variable 
to be still well defined even if z = 0, which is not the case 
for w and 51 individually. Finally, dividing by y^l — e^) 
avoids to test non-physical values at large eccentricities. 

The following equation is then used: 



(1 



sin - sin(w — 51) 



(A.l) 



where vq is the true anomaly of the planet at a specific 
better constraining date to. Here we fix to to be the date 
of the 2003 observation (Nov. 13, 2003). We run 10 chains 
in parallel and use the Gelman-Rubin statistics as conver- 
gence criterion (see details in Ford 2006). The results of the 
MCMC runs are reported in the context of /3 Pic b in Fig. 3 
and 4. 



Appendix B: Orbital fit material 

You will find below the astrometric data of (3 Pic b(with 
their error bars from Table 2), together with the results of 
the orbital fit for: 1/ the best LSLM solution, and 2/ the 
" highly probable orbit" according to the MCMC approach. 
The details of the orbital parameters of both solutions are 
given in Table 3). Fig. B.l gives the results of the orbital 
fit on the projected sky plane. Fig. B.2 and B.3 give the 
evolution of the relative astrometric values of Aa and AS 
as a function of time. 
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Fig. B.l. Measured astrometric positions of /3 Pic b relative 
from A on the plane of sky used in the present work to char- 
acterize the orbital parameters of the (5 Pic b planet. Together 
with the observed measurements are overplotted the orbital so- 
lution the best LSLM model and the highly probable orbit 
obtained with the MCMC approach. 
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Fig. B.2. Measured astrometric positions of j3 Pic b relative 
from A in Aa as a function of time. Both orbital solutions of 
the best LSLM model and the highly probable orbit obtained 
with the MCMC approach are overplotted. 
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Fig. B.3. Measured astrometric positions of j3 Pic b relative 
from A in AS as a function of time. Both orbital solutions of the 
best LSLM model and the highly probable orbit obtained 
with the MCMC approach are overplotted. 



